

% who are you?
% analyst = 1,...,4;

 
%%
load(['sim_eq_round2_',num2str(analyst),'a_replication.mat'])
rng(srng)

% strategy rofiles (0: abstain; 1: support government; 2: support rebels)
YY=zeros(1,6);
for i5=0:2
    for i4=0:2
        for i3=0:2
            for i2=0:2
                for i1=0:2
                    YY=[YY;[1,i1,i2,i3,i4,i5]]; %#ok<AGROW>
                end
            end
        end
    end
end
clear i1 i2 i3 i4 i5


%% set mean utilities
load('parameter_estimates_direction_round1_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
U0 = payoffs(thetaStar,X,Z,YY,reb,mpow); % mean utilities
clearvars -except analyst N specif thetaStar U0


%% simulated equilibria
N1 = floor(N/4);
NS = 3e+3;

if analyst<4
    
    V_R = U0(:,1,N1*(analyst-1)+(1:N1)) + randn(size(YY,1),NS,N1); % entry payoffs of rebel group
    V_US = U0(:,2,N1*(analyst-1)+(1:N1)) + randn(size(YY,1),NS,N1);
    V_UK = U0(:,3,N1*(analyst-1)+(1:N1)) + randn(size(YY,1),NS,N1);
    V_France = U0(:,4,N1*(analyst-1)+(1:N1)) + randn(size(YY,1),NS,N1);
    V_Russia = U0(:,5,N1*(analyst-1)+(1:N1)) + randn(size(YY,1),NS,N1);
    V_China = U0(:,6,N1*(analyst-1)+(1:N1)) + randn(size(YY,1),NS,N1);

    filen=1:NS;
    intervEq=cell(NS,N1);

    for n=1:N1
        tic
        parfor b=1:NS
            fn=filen(b);
            fname=['game',num2str(analyst),'a',num2str(fn),'.nfg'];
            fileID=fopen(fname,'w');
            fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
            fprintf(fileID,'{ "US" "UK" "France" "Russia" "China" } { 3 3 3 3 3 }\n',1); %#ok<CTPCT>
            formatSpec='%12g ';
            fprintf(fileID,formatSpec,reshape([V_US(2:end,b,n),V_UK(2:end,b,n),V_France(2:end,b,n),V_Russia(2:end,b,n),V_China(2:end,b,n)]',1,5*3^5)); %#ok<*PFBNS>
            fclose(fileID);        
            [~,NE]=system(['gambit-enumpure ',fname]);
            findNE=strfind(NE,'NE');
            EQ=zeros(3^5,length(findNE));
            EV_R=zeros(1,length(findNE));
            for m=1:length(findNE)
                sigma_US0=str2double(NE(findNE(m)+3));
                sigma_US1=str2double(NE(findNE(m)+5));
                sigma_UK0=str2double(NE(findNE(m)+9));
                sigma_UK1=str2double(NE(findNE(m)+11));
                sigma_FR0=str2double(NE(findNE(m)+15));
                sigma_FR1=str2double(NE(findNE(m)+17));
                sigma_RU0=str2double(NE(findNE(m)+21));
                sigma_RU1=str2double(NE(findNE(m)+23));
                sigma_CH0=str2double(NE(findNE(m)+27));
                sigma_CH1=str2double(NE(findNE(m)+29));
                EQ(:,m)=(sigma_US0.^(YY(2:end,2)==0)).*(sigma_US1.^(YY(2:end,2)==1)).*...
                    ((1-sigma_US0-sigma_US1).^(YY(2:end,2)==2)).*...
                    (sigma_UK0.^(YY(2:end,3)==0)).*(sigma_UK1.^(YY(2:end,3)==1)).*...
                    ((1-sigma_UK0-sigma_UK1).^(YY(2:end,3)==2)).*...
                    (sigma_FR0.^(YY(2:end,4)==0)).*(sigma_FR1.^(YY(2:end,4)==1)).*...
                    ((1-sigma_FR0-sigma_FR1).^(YY(2:end,4)==2)).*...
                    (sigma_RU0.^(YY(2:end,5)==0)).*(sigma_RU1.^(YY(2:end,5)==1)).*...
                    ((1-sigma_RU0-sigma_RU1).^(YY(2:end,5)==2)).*...
                    (sigma_CH0.^(YY(2:end,6)==0)).*(sigma_CH1.^(YY(2:end,6)==1)).*...
                    ((1-sigma_CH0-sigma_CH1).^(YY(2:end,6)==2));
                EV_R(m)=V_R(2:end,b,n)'*EQ(:,m);
            end
            intervEq{b,n}=[EV_R;EQ];
        end
        toc
    end

else
    
    V_R= U0(:,1,3*N1+1:N) + randn(size(YY,1),NS,N-3*N1); % entry payoffs of rebel group
    V_US = U0(:,2,3*N1+1:N) + randn(size(YY,1),NS,N-3*N1);
    V_UK = U0(:,3,3*N1+1:N) + randn(size(YY,1),NS,N-3*N1);
    V_France = U0(:,4,3*N1+1:N) + randn(size(YY,1),NS,N-3*N1);
    V_Russia= U0(:,5,3*N1+1:N) + randn(size(YY,1),NS,N-3*N1);
    V_China= U0(:,6,3*N1+1:N) + randn(size(YY,1),NS,N-3*N1);

    filen=1:NS;
    intervEq=cell(NS,N-3*N1);

    for n=1:N-3*N1
        tic
        parfor b=1:NS
            fn=filen(b);
            fname=['game',num2str(analyst),'a',num2str(fn),'.nfg'];
            fileID=fopen(fname,'w');
            fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
            fprintf(fileID,'{ "US" "UK" "France" "Russia" "China" } { 3 3 3 3 3 }\n',1); %#ok<CTPCT>
            formatSpec='%12g ';
            fprintf(fileID,formatSpec,reshape([V_US(2:end,b,n),V_UK(2:end,b,n),V_France(2:end,b,n),V_Russia(2:end,b,n),V_China(2:end,b,n)]',1,5*3^5));
            fclose(fileID);        
            [~,NE]=system(['gambit-enumpure ',fname]);
            findNE=strfind(NE,'NE');
            EQ=zeros(3^5,length(findNE));
            EV_R=zeros(1,length(findNE));
            for m=1:length(findNE)
                sigma_US0=str2double(NE(findNE(m)+3));
                sigma_US1=str2double(NE(findNE(m)+5));
                sigma_UK0=str2double(NE(findNE(m)+9));
                sigma_UK1=str2double(NE(findNE(m)+11));
                sigma_FR0=str2double(NE(findNE(m)+15));
                sigma_FR1=str2double(NE(findNE(m)+17));
                sigma_RU0=str2double(NE(findNE(m)+21));
                sigma_RU1=str2double(NE(findNE(m)+23));
                sigma_CH0=str2double(NE(findNE(m)+27));
                sigma_CH1=str2double(NE(findNE(m)+29));
                EQ(:,m)=(sigma_US0.^(YY(2:end,2)==0)).*(sigma_US1.^(YY(2:end,2)==1)).*...
                    ((1-sigma_US0-sigma_US1).^(YY(2:end,2)==2)).*...
                    (sigma_UK0.^(YY(2:end,3)==0)).*(sigma_UK1.^(YY(2:end,3)==1)).*...
                    ((1-sigma_UK0-sigma_UK1).^(YY(2:end,3)==2)).*...
                    (sigma_FR0.^(YY(2:end,4)==0)).*(sigma_FR1.^(YY(2:end,4)==1)).*...
                    ((1-sigma_FR0-sigma_FR1).^(YY(2:end,4)==2)).*...
                    (sigma_RU0.^(YY(2:end,5)==0)).*(sigma_RU1.^(YY(2:end,5)==1)).*...
                    ((1-sigma_RU0-sigma_RU1).^(YY(2:end,5)==2)).*...
                    (sigma_CH0.^(YY(2:end,6)==0)).*(sigma_CH1.^(YY(2:end,6)==1)).*...
                    ((1-sigma_CH0-sigma_CH1).^(YY(2:end,6)==2));
                EV_R(m)=V_R(2:end,b,n)'*EQ(:,m);
            end
            intervEq{b,n}=[EV_R;EQ];
        end
        toc
    end

end


%%
clear filen n
save(['sim_eq_direction_round2_',num2str(analyst),'a.mat'])


